Dominik Adamczyk

Mownit Lab9 rozwiązania

1. Zamień sygnał na sumę sygnałów (1 pkt) np

x = sin.(2*pi*t*200) + 2* sin.(2*pi*t*400)

Zaobserwuj wynik transformaty i wyjaśnij go.

using Plots, FFTW
Fs = 1024
t = 0:1/(Fs-1):1
x = sin.(2*pi*t*200) + 2* sin.(2*pi*t*400)

plot(t, x)

y = fft(x)
sticks((abs.(fft(x))))

Widzimy dwa sygnały dla częstotliwości 200 i 400 zgodnie z danymi w zadaniu. Dodatkowe widoczne sygnały wynikają z natury algorytmu FFT. FFT działa na zbiorze liczb zespolonych, i powiązuje sygnał z zespolonymi sinusoidami. Zespolono sinusoida jest spiralą mającą dwie części - realną i zespoloną. Części te są przesunięte względem siebie o $\pi / 2$. Przez to są ortogonalne i mogą "złapać" komponent sygnału na tej częstotliwości. Sygnały te powstają dla sinusoidy przesuniętej o 9o stopni

https://dsp.stackexchange.com/questions/4825/why-is-the-fft-mirrored#:~:text=Because%20both%20the%20positive%20and,signals%20in%20the%20same%20way.

2. Usuwanie szumów (1 pkt) :

   1. Wypełniamy tablicę wartościami funkcji cosinus ("sygnału") zaburzonej niewielkim "szumem" np. dodając do każdej wartości wylosowaną liczbę funkcją rand().
   2. Proszę narysować wykres zaszumionej funkcji.
   3. Narysować wykres transformaty Fouriera (widmo) tego sygnału (jak poprzednio).
   4. Po transformacie wyzerowac w widmie wszystkie elementy, których wartość bezwzględna jest mniejsza niz 50. W ten sposób usuwamy "szumy" z sygnału.

   5.Przeprowadzić odwrotną transformatę funkcją ifft(). Narysować wykres otrzymanej funkcji dla częsci rzeczywistej. Porównać z wejściowym wykresem sygnału.

Fs = 1024
t = 0:1/(Fs-1):1
noised = cos.(2*pi*t*10) .+ rand(length(t))
plot(t, noised)

y = fft(noised)
sticks((abs.(y)))

filtered = copy(y)

filtered = map(x -> if(abs(x) < 50) 0 else x end, filtered)

plot(t, real(ifft(real(filtered))))

plot(t, noised, label="noised")
plot!(t, real(ifft(real(filtered))), linewidth=2, label="cleared")

3. Proszę nagrać własny glos i zastosować na nim trasformatę Fouriera, narysować wykres widma. Następnie poeksperymentować (wyciąć wybrane częstotliwości), dokonać odwrotnej transformaty i odsłuchać efekt (3 pkt) .

using WAV

s, freq = wavread("voice.wav")
s1 = s[:, 1]
freq = Int(freq)

timeArr = (0:(size(s)[1]-1)) / freq * 1000
plot(timeArr, s1)

f = fft(s1)
sticks(abs.(f))

filtered = copy(f)
for i in 1:length(filtered)
    if abs(filtered[i]) < 20
        filtered[i] = 0
    end
end

plot(real(ifft(filtered)))

wavwrite([i.re for i in ifft(filtered)], "out1.wav", Fs=freq)

filtered = copy(f)
for i in 1:length(filtered)
    if abs(filtered[i]) > 20
        filtered[i] = 0
    end
end

plot(real(ifft(filtered)))

wavwrite([i.re for i in ifft(filtered)], "out2.wav", Fs=freq)

Zaszumianie i odszumianie nagrania

s, freq = wavread("voice.wav")
s1 = s[:, 1]
freq = Int(freq)

s1 .+= (rand(length(s1)) .-1 ) / 10
# wavplay(s1, freq)
wavwrite(s1, "noisyAudio.wav", Fs=freq)

plot(s1)

transformata = fft(s1)
sticks(abs.(transformata))

filtered = copy(f)
for i in 1:length(filtered)
    if abs(filtered[i]) < 13 || abs(filtered[i]) > 80
        filtered[i] = 0
    end
end

wavwrite([i.re for i in ifft(filtered)], "clearedAudio.wav", Fs=freq)

plot(real(ifft(filtered)))